/**
 * Copyright 2012 Brigham Young University
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *    http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package edu.byu.nlp.cluster.mom;

import java.util.logging.Logger;

import org.apache.commons.math3.special.Gamma;

import edu.byu.nlp.cluster.em.ExpectationMaximization.ParameterPair;
import edu.byu.nlp.cluster.em.HyperParameterOptimizer;
import edu.byu.nlp.math.optimize.ConvergenceChecker;
import edu.byu.nlp.math.optimize.ConvergenceCheckers;
import edu.byu.nlp.util.Matrices;

/**
 * @author rah67
 *
 */
public class VariationalBayesHyperparameterOptimizer implements HyperParameterOptimizer<HyperParams, MoMParameters> {

	private final static Logger logger = Logger.getLogger(VariationalBayesHyperparameterOptimizer.class.getName());
	
	/** {@inheritDoc} */
	@Override
	public HyperParams optimizeHyperParameters(MoMParameters curParams, ParameterPair<HyperParams> hyperParamPair) {
		maximizeBeta(curParams, hyperParamPair);
		return hyperParamPair.getNextParams();
	}

	/**
	 * x_{n+1} = x_n - f'(x_n) / f"(x_n)
	 * f'(beta) = KV * digamma(beta * V) - KV * digamma(beta) + \sum_k (digamma( b_kv) - digamma(b_k0))
	 * f"(beta) = KV^2 * trigamma(beta * V) - KV * trigramma(beta)
	 * 
	 * @return
	 */
	private void maximizeBeta(MoMParameters curParams, ParameterPair<HyperParams> hyperParamPair) {
		final int K = curParams.getNumClasses();
		final int V = curParams.getNumFeatures();
		// a and b have now been altered to contains differences of logarithms of digammas.
		double[][] digammaBMinusDigammaB0 = curParams.getLogPOfXGivenY();
		double beta = hyperParamPair.getCurParams().getBeta();
		int i = 0;
		double prev = beta;
		ConvergenceChecker cc = ConvergenceCheckers.or(ConvergenceCheckers.relativePercentChange(1e-6), ConvergenceCheckers.maxIterations(250));
		do {
			double first = K * V * Gamma.digamma(beta * V) - K * V * Gamma.digamma(beta) 
					+ Matrices.sum(digammaBMinusDigammaB0);
			double second = K * V * V * Gamma.trigamma(beta * V) - K * V * Gamma.trigamma(beta);
			double nextBeta = beta - first / second;
			prev = beta;
			beta = nextBeta;
		} while (!cc.isConverged(++i, prev, beta));
		logger.info("Beta = " + beta);
		hyperParamPair.getNextParams().setBeta(beta);
	}

}
